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We use Bayesian statistics and Markov chain Monte Carlo (MCMC) tech- 
niques to construct dynamical models for the spiral galaxy NGC 6503. The 
constraints include surface brightness profiles which display a Freeman Type II 
structure; H I and ionized gas rotation curves; the stellar rotation, which is nearly 
coincident with the ionized gas curve; and the line of sight stellar dispersion, 
which displays a cr— drop at the centre. The galaxy models consist of a Sersic 
bulge, an exponential disc with an optional inner truncation and a cosmologically 

, motivated dark halo. The Bayesian/MCMC technique yields the joint posterior 

probability distribution function for the input parameters, allowing constraints 
>• ■ on model parameters such as the halo cusp strength, structural parameters for 

the disc and bulge, and mass-to-light ratios. We examine several interpretations 
of the data: the Type II surface brightness profile may be due to dust extinction, 
to an inner truncated disc or to a ring of bright stars; and we test separate fits 

q ■ to the gas and stellar rotation curves to determine if the gas traces the gravi- 



tational potential. We test each of these scenarios for bar stability, ruling out 
dust extinction. We also find that the gas likely does not trace the gravitational 
potential, since the predicted stellar rotation curve, which includes asymmetric 
drift, is then inconsistent with the observed stellar rotation curve. The disc is 
well fit by an inner-truncated profile, but the possibility of ring formation by a 
bar to reproduce the Type II profile is also a realistic model. We further find 
that the halo must have a cuspy profile with 7 > 1; the bulge has a lower M/L 
than the disc, suggesting a star forming component in the centre of the galaxy; 
and the bulge, as expected for this late type galaxy, has a low Sersic index with 
n& ~ 1 — 2, suggesting a formation history dominated by secular evolution. 
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Introduction 



The joint analysis of surface brightness profiles and rotation curves has traditionally 

yielded rich insights into the physics of disc galaxies. Surface brightness profi les allow for 

bulge-disc decomposition, while extended rotation curves constrain mass models (Ivan Albada et al 
198,4 IPalunas fc Williamsl l2~00~ol : Ide Blok et allhoOli kuzio de Narav et al.lhoosh . However, 
mass models are subject to degeneracies due, in part, to uncertainties in the mass-to-light 
(M/L) ratios of the baryonic components (e.g., Mailer et al. 2000, Dutton et al. 2005). 
Breaking these degeneracies requires more data, such as a line of sight (LOS) velocity dis- 
persion profile o r a colour profile. Additional data of this type allows for t he creation of dy- 
nami c al models f|Rix et al 1119971 : iGebhardt et al.ll2000l : IWidrow et al.ll2003l : iBaes fc Dejonghe 



2004; Thomas et al. 



20071 ). Dynamical models, in turn, provide the initial conditions for 



N— body simulations, which are needed to study the growth of bar and spiral instabili- 
ties. In addition, non-circular motions co mplicate the interpretation of gas rotation curves 
(jRhee et al.l I2004J : IValenzuela et al.l 120071 ) . Stellar rotation curve data may therefore help 
address these complications. A data set that includes surface brightness profiles, gas and 
stellar rotation curves and stellar velocity dispersions can provide constraints on fundamental 
galaxy parameters, along with information about stability to bar formation in the disc. 

Such a data set exists for the isolated dwarf Sc galaxy NGC 6 503. Gas and s tellar 



rotation c urves were measured by Ide Vaucouleurs fc Cauletl (119821 ) . iBegemanl (119871 ) and 
Bottemal (Il989l ) (hereafter B89); B89 also measured the stellar velocity dispersion profile 
and surface brightness profiles in the B— and R— bands. There are no strong asymmetries 
in galaxy structure (such as a bar), so it is a good candidate for analytic axisymmetric 
models. However, further investigation reveals several peculiarities. The galaxy possesses 
a sharp drop in velocity dispersion near the centre (known as a 'a— drop'). The gas and 
stellar rotation curves are nearly coincident despite the standard asymmetric drift formula's 
prediction that stars should noticeably lag behind the gas. Finally, the surface brightness 
profile displays four distinct regions: a central peak, a flat region indicating a Freeman 
(1970) Type II profile, an inner exponential of scale length ~ 1 kpc and a shallower outer 
exponential. 

In this paper, we present several scenarios corresponding to differing interpretations of 
the data. For each scenario, we construct dynamical models that reproduce these features. 
We decompose the surface brightness using a Sersic bulge and an inner-truncated light profile 
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of the type used by iKormendyi (Il977l ). Inner truncated disc p rofiles may be due to dust 
extinction against an exponential disc ( iMacArthur et al.l 120031 ) ; they may be intrinsic to 
the density profile; or they may be due to a ring of bright stars which does not trace the 
mass. We explore these possibilities by testing models numerically under each scenario for 
bar formation. The coincidence of gas and stellar rotation curves also deserves attention as 
it suggests far less asymmetric drift than predicted. Thus, we wish to understand whether 
gas follows circular orbits that trace the gravitational potential, or whether it has its own 
asymmetric drift. We investigate this issue by testing two different ways of fitting the rotation 
curves: by fitting the stellar rotation alone and using the asymmetric drift to calculate the 
circular velocity; and by fitting the gas alone assuming it traces the circular velocity. 

Our aim is to construct dynamical models for each of these scenarios and test them via 
N— body simulations. We use the GalactlCS model from Widrow et al. (2008, hereafter 
WPD) defined in terms of distribution functions for the disc, bulge and halo. The halo 
allows for a cusp or a core, and the bulge follows a Sersic (1968) profile. The disc surface 
density profile is exponential with an optional inner truncation. 

The large number of parameters makes finding a fit to the data a challenging exercise. A 
number of techniques have been developed to determine best-fit parameters in such complex 
spaces - notably, maximum likelihood techniques which involve minimizing the y 2 value. 
One promising approach which improves upon the maximum likelihood method employs 
Bayes' theorem and a Markov chain Monte Carlo (MCMC) technique to survey the relevant 
parameter space. Bayes' theorem provides a method of determining the posterior probability 
of a hypothesis based on both the likelihood function and prior information about the input 
parameters; MCMC provides for a way to survey the parameter space. In conjunction, the 
two techniques provide the joint posterior probability distribution function (PDF) of the mul- 
tidimensional parameter space; the marginal posterior PDF for any parameter (and hence its 
formal mean and error bars) are obtained by marginalizing over nuisance parameters. This 
technique yields realistic constraints on the model parameters, including (but not limited 
to) the halo cusp, the disc and bulge M/L ratios, masses for each component, inner trun- 
cation parameters, and bulge Sersic ind ex. Similar techniques have been used, for example, 
to de termine cosmological parameters (ITegmark et al.l 120041 ; iPercivall 120051 ; ICorless fc King 
2008T ). 



In this work we build dynamical models for NGC 6503 using a Bayesian/MCMC method. 
We fit several data sets for this galaxy simultaneously in order to obtain as complete a picture 
of the galaxy as possible; however, the lack of a bar provides a further constraint. There- 
fore, we also obtain the stability parameters X a nd Q for our m odels. The X— parameter 
determines global stability to multi-armed modes (lToomrelll98ll ) and the Q— parameter de- 
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termines stability locally ( lToomrdll964r). As a purely phenomenol ogical matter, Q can also 
determine global stability to bars ( lAthanassoula fc Sellwoodlll986l ). We test the stability of 
our models using N— body simulations. 

Our approach allows us to investigate numer ous other properties of the galaxy. The 
Sersi c index is indicative of bulge formation history ( Courteau et al.|[l996 : Kormendy fc Kennicutt 
20041 1 . We also explore the possibility that the bulge is related to luminous nuclear clusters 
found near the centres of many disc-dominated late-type galaxies, and derive Mj L ratios for 
the disc and bulge. Moreover, we can constrain the cuspiness of the halo. Cosmological sim- 
ulations consistently demonstrate that halos are expected to have cuspy halos (e.g., Navarro 
et al. 1997; Moore et al. 1999), but there are considerable difficulties involved in inferring 



central halo densities from rotation curve mea surements (IHayashi et al.l I2004J : iRhee et al. 



20041 ; iDutton et al.ll2005l ; IValenzuela et al.ll2007l ). For example, gas rotation curve data near 
the centre of the galaxy may not trace the gravitational potential; using th e stellar rota- 



tiona l velocity, where available, may be a better way to infer the cusp value (IPizzella et al. 



20081 ). Further, noncircular gas motions at the centre (as may be caused by a triaxial halo) 
may skew the inter pretation of the rotation curve, depending on the orientation of the disc 
(jSimon et al.l 120051 ). With both gas and stellar rotation curves, NGC 6503 provides an 
excellent opportunity to investigate these issues. 

This paper is organized as follows. In §2 we detail the observational properties of NGC 
6503 and elaborate on the peculiarities outlined above. In §3 we discuss our dynamical 
models, and in §4 we describe our application of the Bayes/MCMC method. We discuss 
our primary stability results in §5, in the process constraining the available physics, and we 
present additional results in §6. Section 7 is the discussion, which includes an exposition of 
prior research on NGC 6503. 



2. NGC 6503: Observations and Issues 



NGC 6503 is an isolated dw arf Sc galaxy at a distance of 5.2 Mpc (IKaranchentsev fc Sharina 



19971 ) and is inclined at 74 ± 1° (lBegemanlll987l ). It displays a small, bright central bulge, no 
signs of strong asymmetric central structure (e.g., no obvious bar) apart from some spirality, 
and it is mostly free of gas. In this section we review existing photometric and kinematic ob- 
servations of the galaxy. Throughout this paper, we use a distance at 5.2 Mpc and therefore 
1 kpc = 40" at this distance. 
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2.1. The surface brightness profile 



Photometric imaging by B89 provides azimuthally averaged surface brightness (SB) 
profiles in the B— and R— bands, both of which display a sharp rise within 50 pc and a 
Freeman (1970) Type II hump between 50 pc and ~ 2 kpc. The SB profile is exponential 
beyond that. The central flattening is accompanied by a slight reddening in the B — R profile 
(Fig. CD) . A circumnuclear Hct ring is observed at R ~ 1.2 kpc, indicating a star-forming 
region (IKnapen et al.ll2006l ). The outer exponential is remarkably straight in spite of the 
large error bars; a linear fit to those points yields a photometric scale length R^ of 1.18 ±0.04 
kpc, which becomes about 1.3 kpc when corrected for inclination and scale height. Fig. [2] 
shows the log //-log R plot of the surface brightness. The inner few points are almost straight, 
approaching a line with a slope ~ 1, suggesting a Sersic index of 1. 

There are three possibilities for the origin of Type II profiles: first, they may be caused 
by dust extinction; second, they may reflect the formation of bright stars that do not trace 
the mass; and third, they may be intrinsic to the mass distribution of the disc. In NGC 6503, 
the reddening in the inner portion may support the dust extinction hypothesis (Bottema & 
Gerritsen 1997, hereafter BG97). BG97 adopt an exponential disc and assume that the 
flattening in the surface brightness is due to dust (we return to this wor k in §7.1). The 
possib ility that dust extinction causes Type II profiles is also discussed by MacArthur et al. 
(120031 ) for their sample of spiral galaxies with optical and IR imaging; these authors are 
unable to reliably distinguish between Type I and Type II profiles on the basis of dust 
extinction, however. 



The second possibility relates to ba r formation. iBaggett. Baggett fc Anderson! (119981 ) 
and lAnderson. Baggett &: Baggettl (120041 ) show that many galaxies with Type II profiles can 
be fit using Kormendy's (1977) inner truncated profile 



I(R) = I ex P [-(R/R d + (R h /R) a )} 



(1) 



where Rh is a turnover radius and a is a cutoff index. Most of Baggett et al.'s (1998) 
galaxies have a bar, but a significant minority do not. These authors also argue that such 
flattening is intrinsic to the galaxy rather than being caused by dust. It is known that 
bars can induce ring formation in galaxies, which may contribute to the Type I I profile 
flatte ning;, and that circumnucl ear rings correlate strongly with barred galaxies (IKnapen 
20051 ). Moreover, simulations by lFoyle et al.l ( 120081 ) suggest that the SB flattening could also 
be caused by the presence of a central bar, as the resulting reorganization of light can flatten 
the SB profile over the bar's extent. Furthermore, while NGC 6503 does not appear to have 
a bar, it is possible for bars to dissolve by gas inflow, which transfers angular momentum 
from the gas to the bar (IBournaud et al.l 120051 ) . The ring currently seen in NGC 6503 may 
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therefore be a relic of a now dissolved bar. In a scenario of ring formation by bar destruction, 
the outer exponential represents the 'true' disc, and the Type II hump i s mainly due to the 
formation of bright stars that do not trace the mass. The results in iFoyle et al.l (120081 ) 
suggest that the Type II profile could be reproduced by the emergence of a bar in a disc of 
scale length 1.3 kpc. Therefore, it is incumbent to investigate the SB profile by fitting only 
the points exterior to 2 kpc. 

A third possibility is that Type II profiles are caused by a genuine mass deficit in the 
inner region, a nd therefore we also fit the full SB profile directly by assuming an inner 
truncated disc ( jKormendyl 1 19 771 ) . This naturally produces a better fit to the hump than to 
the outer points; along with the first scenario, this scenario provides a direct comparison 
with the modelling of BG97, who did not account for the outer exponential. 



2.2. Rotation curves 



Several rotation curves for NGC 6503 have been measured. H I observations by lBegeman 



(119871 ) reveal a remarkably flat rotation curve between 3 kpc and 20 kpc de viating by no 
more than ~ 4 — 5%, and a maximum circular velocity of ~ 120 km s -1 . iGreisen et al. 



(120091 ) find a similar H I curve in addition to evidence for a thick and thin H I disc in NGC 
6503. The thick disc appears to rotate more slowly than the thin disc and the former does 
not extend beyond the optical radius. However, this distinction does not affect Begeman's 
measurements of the total H I rotation curve. B89 maps the inner rotation curve w i th H/3 
and O III emission as well as the stellar features, while Ide Vaucouleurs fc Cauletl (119821 ) 
supply Ha data. Where they overlap, all data sets are in agreement; the H I observations 
appear to be slightly but systematically larger than the H/3. The stellar rotation curve 
largely coincides with the H/3 and O III curves, with an average difference of only ~ 2 km 
s _1 . At about 1.5 kpc, the stellar rotation exceeds the gas rotation, possibly owing to the 
effect of a spiral arm. For our models, we fit the H I, H/3 a nd stellar data sets. These are 
found in Fig. [31 along with fits based on the fitting formula (ICourteaulll997l ): 



v(r) 



v 



(2) 



[l + (R c /R)<y/< 

where vq is an asymptotic velocity, R c is a projected length scale and ( is a shape parameter 
governing the sharpness of the 'turnover' to the asymptotic velocity. The fit suggests a 
projected asymptotic circular velocity of 117±1 km s _1 and shows how close the two rotation 



curves are. 



The coincidence of the stellar and gas rotation curves must be accounted for. If the gas 
rotation traces the gravitational potential, the difference between the two rotation curves 
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yields a measure of the asymmetric drift, v a . The latter, defined as the difference between 
the circular velocity v c calculated from the potential and the streaming velocity v s of the 
stars, is given by 



a 



Va 



R 

2Vr 



2R 1 
Rd 2 



R dv s 
v s dR 



(3) 



for an exp onential disc, where Rd is the disc scale length and ctr is the radial velocity 
dispersion (jBinney fc Tremaind 120081 ). When applied to the observed H/3 rotation curve, 
Eq. [3] yields an asymmetric drift of order of 5 — 15 km s _1 , depending on the adopted 
central radial dispersion^] The asymmetric drift is much larger than the observed difference 
between the gas and stellar curves; a possible cause of the discrepancy is that gas does not 
trace the gravitational potential. Traditionally, gas is assumed to travel on circular orbits, so 
that the H/5/H I rotation curves would trace the gravitational potential in this galaxy. The 
assumption of circularity in gas rotation curves underlies much of the work on mass modelling 
in the literature. However, it is not possible for this assumption to strictly hold because gas 
dispersion is non-zero. There are several sources of this dispersio n including turbulenc e 
caused by supernova feedback and the magnetorotational instability f lTamburro et al.ll2009l ). 
In addition, non-circular motions can impact the gas rotation curves of disc galaxies (e.g., 
Hayashi et al. 2004, Valenzuela et al. 2007), thereby reducing the observed velocity relative 
to the actual circular velocity. Being collisionless, our models cannot capture these effects; 
therefore, the hypothesis that the gas rotation traces the circular velocity can only be tested 
using the gravitational potential of our models. However, the stellar asymmetric drift can 
easily be modeled to obtain a separate value for the circular velocity. The stellar asymmetric 
drift is itself subject to several caveats; it is, for example, highly sensitive to the orientation 
of the velocity ellipsoid. These issues are addressed in §7.3. 



2.3. The LOS velocity dispersion 

The LOS velocity dispersion profile for NGC 6503 from B89 is found in Fig. HJ There 
is a cr-drop within 300 pc and an exponential decline beyond that range. We fit the data to 
an outer exponential with an inverse Gaussian in the centre: 

a{R) = a exp(-R/R a ) [l - Aexp(-R 2 /(2B 2 ))] (4) 



1 T589 claims that the calculated asymmetric drift implied by the rotation curves is ~ 3 — 5 km s -1 , which 
is much lower than our calculated values. We have been unable to reproduce this number; every permutation 
of scalelength and dispersion that we tried in Eq. [3] yields asymmetric drifts of at least ^5 km s _1 . 
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where A, B, the extrapolated central dispersion a Q and the projected scale length R a are the 
fitting parameters. The fit yields R a = 1.16 ± 0.29 kpc and cxq = 59.5 ± 8.4 km s , which 
will provide a point of comparison for our results in §6. Note that we do not attempt to 
model the physical processes which may be responsible for the a— drop, but we do determine 
structural parameters and effective M/L ratios for the bulge to account for the a— drop. 

NGC 6503 was the first galaxy shown to have a a— drop. A growing body of data suggests 
that a— drops are found in numerou s spiral g alaxies (|Marquez et al. 2003 : Comeron et al. 



2008 ; de Lorenzo-Caceres et alJ[2008 ). although Kovela et al. ( 20081 ) argue that some a— drops 



are observational artifacts. BG97 suggest that the cr— drop is due to a small central, possibly 
star forming component ca pable of continually c oolin g the centre. A physical explan ation 
for a— drops is proposed by IWozniak et al.l (120031 ) and lWozniak fc Champavertl (120061 ) . who 
argue that the presence of dynamically cold gas funnelling into the centre is responsible for 



formin g dynamically cold stars that dominate the observed LOS kinematics. IComeron et al. 



( 120081 ) 's sample of a— drop galaxies correlates with higher incidences of nuclear dust spirals, 
Hoj rings and Seyfert fraction, which would be co nsistent with a cold star-f o rming compo- 
nent in the centre due to gas inflow. Meanwhile, Ide Lorenzo-Caceres et al.l (120081 ) suggest 
that inner nuclear bars may be responsible for at least some of the a— drops found in the 
literature. 



2.4. Four scenarios to test 

We test three different ways to fit the SB profile: (i) assuming an inner truncated disc; 
(ii) assuming an exponential disc with an inner truncated light profile to model simply the 
effect of dust extinction (a full dust model is outside the scope of this work) ; and (iii) fitting 
only the inner [R < 0.1 kpc) and outer {R > 2 kpc) points to account for the true disc being 
external to the Type II hump. Furthermore, we test two ways to model the rotation: (i) we 
fit the gas rotation by assuming that it equals the circular velocity; and (ii) we fit the stellar 
rotation by simulating the LOS observations of the inclined disc. Simulating observations 
along a line of sight is particularly important for stellar kinematic measurements because of 
the non-negligible scale height of the stellar disc, which leads to integration effects as the 
line of sight passes through regions of the disc that are not in the plane. This effect is less 
important for the gas disc, which is expected to have a small scale height. 

These tests are accounted for by four different scenarios. In scenario K (for Kormendy), 
we fit the full surface brightness with an inner truncated Kormendy disc, and we fit the 
stellar rotation curve. In KL (Kormendy light), we fit the surface brightness using an inner 
truncated light profile atop an exponential disc to account for dust extinction, and fit the 
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stellar rotation as in scenario K. In scenario KG (Kormendy gas), we fit the surface brightness 
as in K but fit the gas rotation by assuming that the gas travels on circular orbits, ignoring 
the stellar rotation data. Finally, in scenario E (exponential), we fit the surface brightness 
excluding the points that constitute the Type II hump, and fit the stellar kinematic data as 
in K. Scenario E therefore alleviates the need for an inner truncated disc. In all cases, we fit 
the LOS velocity dispersion and the H I data outside 3 kpc. These models and the fitting 
procedures are described in more detail in §4 and 5; we first briefly discuss the GalactlCS 
model and the MCMC technique. 



3. GalactlCS Model 



Our goal is to construct a self-consistent equilibrium dynamical model for NGC 6503. 
Dynamical modelling yields the phase space distribution function (DF) for a system that 
precisely specifies its densities and velocities simultaneously; for self-consistency and equilib- 
rium, the model must satisfy both the Poisson and time-independent collisionless Boltzmann 
equations. The GalactlCS model of Widro w et al. (2008, hereafter W PP) is designed to 



satisfy these criteria; it is derived from the iKuijken fc Dubinskil (11995) model. We briefly 



describe the GalactlCS model here and refer the interested reader to 



Kuijken fc Dubinski 



dl995|) and WPD for details. 



3.1. Model components 



3.1.1. The bulge and halo 



The bulge is designed to reproduce the Sersic profile (ISersid Il968l ; ICiottil Il99ll ) upon 
projection, for which the intensity is given by 



J(r) = J exp [-b(r /r b ) 1/n »] 



(5) 



where r is the projected radius, J is the central intensity, r b is the half-light radius, n b is 
the Sersic index of the bulge and b is a parameter which depends on n b . The expression is 
deprojected using an Abel integral equation to yield the intrinsic density distribution 



Pbulge(r) 



r 

Pb I — 

Jb 



exp 



-b\ r - 

n 



l/n 



(6) 
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where r is the spherical r adius and p — 1 — 0.6097/ra + 0.05563/n 2 yields the Sersic profile 
( iPrugniel fc Simienl 119971 ). We define a characteristic velocity scale Vb for the bulge by 



[Aixnb n ^T(n(2 - p))r 2 bPb \ 



1/2 



(7) 



and use this as an input parameter instead of pb- 



For the halo we adopt a cosmologically motivated profile (e.g., Navarro, Frenk & White 
1996, 1997; Moore et al. 1999; Diemand et al. 2005), given by 



P 



(r/a fc )T(l + r/a h f-< 



where ph and are the halo scale density and radius, respectively, and 7 is the cusp value 
governing the shape of the halo; 7 = 1 yields a NFW profile, while 7 = produces a core. 
A characteristic velocity scale 0^ is introduced according to o\ = Ana^ph. For practical 
purposes, the density profile is truncated using an error function (WPD). We allow 7 to 
vary for maximum flexibility in fitting the observed data, si nce simulations of halo formation 
suggest that cusp values other than 1 may be appropriate (|Moore et al.lll999l) and that halo 
profiles may not be universal (jJing k, Sutoll2000t iNavarro et al.ll2008l ). 



3.1.2. The disc and the Kormendy profile 

Surface brightness profiles in discs tend to decline exponentially with radius. Assuming 
light traces mass, this implies that the surface density is exponential: 



£(#) = S exp(- J R/i? fl 



(9) 



where Rd is the photometric scale length. Additionally, observations of edge -on galaxies sug- 
gest t hat disc galaxies have vertical sech 2 profiles with constant scale heights ( Ivan der Kruit &: Searle 
198ll ). Thus, the three dimensional density of the GalactlCS disc is given by 



p{R,z) = p exp(-R/R d )sech 2 (z/z d ) 



(10) 



where z& is the scale height. The DF corresponding to this density is given in lKuijken &: Dubinski 
( 119951 ) : it is an extension of Shu's (1969) planar DF into three dimensions. 



However, as discussed in §2, Type II profiles can also be modelled by Eq. [TJ If the light 
traces the mass, this suggests a surface density given by 



E{R) = S ex P [-{R/R d + {R h /R) a )\ . 



'IF 
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As before, the three dimensional density of the disc is given by 



p(R,z) = Po exp [-{R/R d - {R h /R) a )]sech 2 {z/z d ). (12) 

We modify the GalactlCS model accordingly to generate this density. If Type II profiles are 
instead assumed to be due to dust extinction, Eq. [10] holds for the disc density, and only the 
light follows an inner truncated profile (Eq. [Q. 



3.1.3. The dynamics 



The KD disc employs two integrals of motion: the energy E and the angular momentum 
L z , and an approximate third integral describing the energy in the vertical motions, E z . By 
the Jeans theorem, any f unction of three isolating i ntegrals of motion in a given potential will 
exactly solve the CBE (IBinney &: Tremaindl2008l ). The third integra l allows for <j z ^ er^; 
observations of the solar neighbourhood suggest that cr z /aR = 0.6 (jWielenl 1 19741 ). which 
would be impossible if the DF were a function of only two integrals of motion. 

The GalactlCS model decouples the vertical and radial dispersion of the disc. The 
vertical dispersion is given by the vertical potential gradient and scale height. The radial 
velocity dispersion is given by 

a 2 R = a 2 exp(-R/R a ) (13) 
and the tangential dispersion is given by the epicycle equations, i.e. 



(14) 



Here, Q = v/R is the angular rotation speed and the epicycle frequency k is given by k 2 = 
[R(dfl 2 /dR) + 4Q 2 }. By design, the velocity ellipsoid in the GalactlCS disc is cylindrically 
aligned. 



Note that Eq. [13] includes a parame ter to control th e scale length of the velocity dis- 
persion profile. Although observations by iBottemal (119931 ) suggest that R d = R a) there is no 
clear theoretical reason why this should necessarily be the case. Here we treat R d and R a 
independently to provide maximum flexibility in fitting all data sets. 



3.2. Summary of the GalactlCS model 

The GalactlCS model, modified to include inner truncation, has two further parameters 
for the disc; the turnover radius Rh and the cutoff index a. We complete the parameter 
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input set with R— band M/L ratios for the bulge and disc. We assume that the M/L ratios 
are constant with radius, and use the same value for the surface brightness and the stellar 
kinematic fits. Although the band pass over which the stellar observations were taken are 
not identical to the R— band, we do not expect the error introduced by this simplification to 
be significant. 

Thus the parameters that we modify are: the five halo parameters; the exponential disc 
mass, scale height, scale length, and two Kormendy parameters; all bulge parameters; the 
central radial dispersion and dispersion profile scale length; and the disc and bulge M/L 
ratios. In addition, each data set has an associated noise parameter that is allowed to vary 
(see §4). We assume the halo to be nonrotating and truncate the disc at 5 kpc (well outside 
the measured data). The galaxy inclination % is also fixed at 74°. A list of the relevant 
parameters that are fit, as well as calculated output quantities, is presented in Table 1. The 
units used by GalactlCS are such that G — 1, but the masses have been converted back to 
10 9 M Q for simplicity. 



4. Statistical Approach 

4.1. Bayes' Theorem and MCMC 

The core idea of this study is to derive the joint posterior probability distribution func- 
tion (PDF) of the galaxy parameters using Bayes' theorem and MCMC. Bayes' theorem 
provides a way to infer the posterior probability of a hypothesis, given some kind of evidence 
(such as observational data), using the prior probabilities of the evidence and model and 
a likelihood function (such as a \ 2 function). The joint posterior probability distribution 
function (PDF) of the model parameters is written as P(M\D), where M is the set of model 
input parameters and D is the set of observational data. Bayes' theorem may be written as 

P( , m = WEM (15) 

where P(M) is the prior probability of the input parameters, P(D\M) is the likelihood of 
the data given a specific set of model parameters and P(D) is the prior probability of the 
data and functions as a normalization. 

MCMC, in turn, provides a method of surveying the parameter space that rapidly 
converges to the posterior probability distribution of the input parameters. Our MCMC 
technique explores the parameter space by way of the Metropolis-Hastings algorithm. Es- 
sentially, a random walk is constructed in parameter space by sampling a new set of model 
parameters M* at each step of the chain, and P(M*\D) is calculated using the likelihood 
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function and the priors. The ratio s = P(M*\D) / P(M\D) is then calculated. If s > 1, M* is 
accepted; otherwise, M* is accepted or rejected with a uniform random probability g. Over 
the course of many such iterations, the chain populates the posterior PDFs, in the process 
supplying the PDFs for each parameter by marginalizing over the full PDF. 

Initially, the chain quickly travels to regions of 'good' fits from any set of initial param- 
eters, and eventually settles into the best fit region, a process known as 'burn- in'. It can be 
proven that the resulting density of points i n parameter sp ace samples the real probability 
density distribution of the parameter space ([Gregory! 120051 ) . We aim towards an acceptance 
rate of 23% in our chains (Roberts et al. 19971 ). A detailed description of the technique may 
be found in I Gregory! (120051 ) . 



The benefits of MCMC are numerous: it fully samples the PDFs for all fitted parameters, 
meaning that formal errors and error contours are easily determined; it is easily extended 
to include more parameters if desired; it tends to move to a region of high likelihood fairly 
rapidly; other quantities, such as the stability parameters Q and X, are easily calculated 
from the PDFs; and there is some e yidence that it pr ovides more realistic error constraints 
than maximum likelihood methods (IKelly et al.l 120081 ) . 



The primary disadvantages are twofold. First, there is no guarantee of uniqueness to 
the final solutions - multimodal PDFs are possible in any given parameter space. However, 
a judicious choice of the step size and proper tweaking to obtain the correct acceptance 
rate mitigates this problem. Second, although MCMC chains tend to converge quickly, it 
is not so simple to verify that the converged chains have fully populated the joint posterior 
PDF. Moreover, the chain may 'spread around' the parameter space, reflecting either the 
real structure of the joint PDF or the fact that the chain has not run long enough to fully 
populate the PDF. Formally ensuring convergence may require prohibitively long computing 
time. However, in practice we are able to obtain good constraints for most parameters in a 
reasonable time. This does not preclude the possibility that a MCMC chain passed through 
multiple PDF peaks after burn-in, but there would be no reason to prefer any one such peak 
over another. 



4.2. The likelihood function and the priors 

The likelihood function is given by 

~(di - mi) 2 ' 



C= (2 
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where N is the total number of data points, d{ is the experimental data value, rrij is the 
mean value predicted by the model and o~i is the error, scale-modified as described below. 
This expression is applied to each data set, and the total likelihood is given by the product 
of the individual likelihoods. 

For the data values cfj, we use the H/3 curve, the LOS dispersion, stellar rotation curve 



and the R— band SB profile from B89 and the extended H I curve from lBegemanl (119871 ). We 
ran several different MCMC chains (described in the next section) making use of different 
data sets. 

The model values are calculated as follows. For the gas, we assume that the observed 
H/3/H I rotation curves follow circular orbits so that the model gas rotation curve can be 
calculated directly from the potential. For the stellar rotation curve and dispersion profile, 
we sample the disc and bulge DFs along the line of sight at each data point from B89, 
appropriately weighted by the Mj L ratios. Finally, for the surface brightness we sample the 
DFs along elliptical annuli around the inclined galaxy centre at an axis ratio determined by 
the inclination angle (i.e., q = cost = 0.28) to determine the elliptically averaged density 
and thus the surface brightness in mag arcsec -2 . We emphasize that in no cases have we 
corrected the observed profiles for inclination; all inclination effects are accounted for by 
rotating the model galaxy and sampling along the line of sight. 

For each dataset used in a particular MCMC run, a single noise parameter is imple- 
mented. This consists of an extra term added in quadrature to each error bar: o\ = + f 2 , 
where is the original experimental error and / is the noise parameter, taken as constant 
for all points in a given data set (so there is a total of four scale factors). This approach 
allows us to account for effects that cannot be captured by our models - for example, that 
of spiral structure on the light profile when the original models are axisymmetric. Noise 
parameters use a Jeffreys prior: 

P(f) = f] ( , 1 /f y (17) 

(Gregory 2005) where f m \ n and / max refer to the minimum and maximum values of the noise 
parameters adopted for the MCMC run. Otherwise, P(M) is uniform for all parameters. 
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5. Results 

5.1. Dust extinction and inner truncation 

As discussed in §2, the flattening in the SB profile may reflect the intrinsic mass dis- 
tribution in the disc, or it may be due to dust extinction. We conduct two MCMC chains 
to examine these possibilities. For the first run, corresponding to scenario K, the mass dis- 
tribution of the disc is given by Eq. [12J so that the disc is intrinsically truncated. We fit 
the full R— band SB profile from B89. In the second run, corresponding to scenario KL, the 
mass distribution of the disc is given by Eq. [10] so that the disc is purely exponential and 
the light follows an inner truncated profile. For both runs, we fit the full stellar rotation 
curve, the H I curve outside 3 kpc and the LOS stellar dispersion profile; we do not fit the 
gas rotation internal to 3 kpc here. 

Fits for representative examples of each scenario are found in Figs. [5H71 In Fig. [5] we 
can see that the models reproduce the stellar rotation well. As expected, the circular velocity 
for these models is always greater than the FL3 data. The rotation curve breakdown is also 
seen here, showing that the inner rotation curve for scenario KL is strongly disc-dominated. 
In addition, the surface brightness fit in Fig. [H] is excellent for both runs within 2 kpc. The 
outer exponential is virtually ignored by the fit owing to the larger error bars. The dispersion 
profile fits in Fig. [7] are also good. Scenario K displays a slightly inferior dispersion profile 
fit, owing to the lower do for these models, but the a— drop is present. The outer region 
fit appears to overshoot the points slightly, owing to the large R a for these runs. We have 
verified that decreasing R a by 50% produces better fits to the dispersion without affecting 
the remaining fits, suggesting that the cause of the poorer dispersion fit in the outer region 
is a local PDF peak. The \ 2 values for scenario KL are approximately equal to those for 
scenario KL; therefore, there is no reason to reject this scenario on the basis of the fit quality. 
However, analysis of the PDFs and the time series trace of the x 2 value reveals that this 
run actually passes through two separate minima which correlate with two distinct regions 
of parameter space, both of which produce nearly identical fits. 

Fig. [8] shows the two-dimensional PDF for ikfj^l and the minimum disc Q } Q m - m , usually 
attained at around 1 — 1.5 disc scale lengths. Scenario K is in black and scenario KL is in 
green. For scenario K, we find that the PDF spans a large range in Md and Q m i n , forming 
a relatively narrow trough from Q m \n ~ 2.3 at Md ~ 2.1 x 1O 9 M in the upper left to 



2 Note that the input M e parameter is the exponential disc mass. It is distinguished from the true disc mass 
Md because the exponential disc must be truncated in the outer region, and because the inner truncation in 
scenarios K and KG removes additional mass from the inner region. 
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Qmin ~ 1-3 at Md ~ 4.4 x 10 9 M Q in the lower right. Not surprisingly, for scenario KL is 
much larger, since there is no hole in the mass distribution. Since the radial dispersion <jr 
and epicycle frequency k are constrained by the data, the increased surface density yields 
lower values for Q mm . For scenario KL, the best fit region extends below the Qmin = 1 line, 
which suggests that some of the best fit models are locally unstable. 

The PDF in Q m i Q and X min is found in Fig. [9l For scenario K, this plot shows that X min 
lies between 1 and 3, and appears to correlate weakly with Q m in, in the sense that models 
with high Q m i n also have high X min . Fig. M also shows that all models from scenario KL 
have X min < 1, suggesting very strong instability to global modes. 

To discern whether or not the models for scenario KL are realistic, we must determine 
whether or not galaxies with Q < 1 are possible. Unstable galaxies generate spiral structure 
that heat stars as they propagate outwards but continual cooling could occur if gas infall 
onto the s urface of the galaxy were substantial enough. Such gas would contribute to star 
formation. iFuchsl (119991 ) conducted a maximum disc analysis for NGC 6503 suggesting that 
Q < 1 implies a star formation rate o f 40 yr -1 ; this cont radicts the rate of 1.5 M Q yr _1 



obtained from the observed Ha flux (IKennicutt et al.lll994l ). Therefore, it is unlikely that 
this cooling method is important in NGC 6503. In the absence of star formation, a disc in 
which Q < 1 may fragment as small scale perturbations grow due to the Toomre instability. 
The likely evolution for such models is to eventually recombine into a single stable disc, but 
not before the stellar dispersions have increased dramatically. 

We select several models from the PDF for both runs (identified by stars on Fig. [8]) 
and evolve them forward in time to exam ine their st ability. To do this, we employ Dehnen's 
(2000) iV-body algorithm implemented by lStifil (120031 ). a fast multipole method that produces 
nearly O(N) scaling. All models are evolved for 5 Gyr with a time-step of 0.5 Myr and a 
softening parameter of 25 pc. We use 500K particles for the disc, 50K for the bulge and 1M 
for the halo. The bar evolution is quantified using the magnitude of the second Fourier mode 
A, which measures the strength of two-armed asymmetries (e.g., Shen & Sellwood 2004): 



.4 



3=1 



N 



(18) 



where 6j is the azimuthal coordinate of the jth particle in the disc. The bar evolution is 
found in Fig. [lOj 

The results show the regions of the M d — Qmin plane that are most susceptible to bar 
formation. For scenario K, we find that much of the region is at least mildly bar unstable, 
with the lower right region more strongly unstable (note that the red curves in Fig. [TUl 
correspond to the discs with lowest Q mm i n each run) and the upper left only showing a mild 
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bar instability. The growth of most bars in this scenario is gradual, with bars only starting 
to become discernable after ~ 2 Gyr. Additionally, the central hole gradually vanishes in 
the bar unstable models as stars stream into the centre during the process of bar formation. 

By contrast, the models for scenario KL are far more susceptible to global instabilities 
than for scenario K, as the growth of the bar mode is rapid and nearly instantaneous. We 
therefore rule out models from scenario KL as being unable to properly model the galaxy; 
this conclusion is a setback for the hypothesis that the Type II SB profile is due to dust 
extinction. 



5.2. Testing the rotation curve fits 

We next test fits to the gas and stellar rotation curves via a third MCMC chain corre- 
sponding to scenario KG. We fit the H/3/H I rotation curves by assuming that they trace the 
circular velocity and ignoring the stellar rotation data. The dispersion and surface bright- 
ness are fit as in scenario K. The stars will rotate more slowly because of asymmetric drift 
(Eq. [3]), an effect that is incorporated into the GalactlCS disc. 

Examining the rotation curve fits in Fig. El it is apparent that scenario K cannot re- 
produce the H/3/H I data and scenario KG cannot fit the stellar rotation data. Thus, it is 
likely that scenario K represents a more realistic model for NGC 6503. Note that both runs 
match the H I data beyond ~ 3 kpc, since there is no asymmetric drift at those radii. The 
other fits, shown by the surface brightness in Fig. [6] and the LOS dispersion in Fig. [7J are 
more than satisfactory. 

Fig. [8] shows that scenario KG displays considerable overlap with scenario K except in 
the higher mass range. The plot reveals that the best fitting models for this scenario span 
a large range in Q m i n , and here some models (albeit very few) extend down to below Q — 1. 
However, Fig. ED shows that all models for scenario KG are confined to X min < 2, unlike for 
scenario K where many models reside between X min = 2 and X min = 3. Also unlike scenario 
K, Qmin does not appear to be significantly correlated with X m i n , and unlike scenario KL, in 
neither scenario K nor KG do we see X m i n drop below 1. 

As before, we selected several models from the PDF for scenario KG and evolved them 
forward in time. The results are found in the bottom panel of Fig. [TTJl As for scenario 
K, most of the available parameter space is bar unstable, with the low mass-high Q part 
of the plot susceptible only to mild bar formation. Thus, it is not possible to distinguish 
between runs K and KG on the basis of bar stability alone. Referring back to the rotation 
curve fits, we see that no models for scenario KG can reproduce the stellar rotation curve 
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because the asymmetric drift is too large. Therefore, we reject models from scenario KG as 
being inconsistent with the data; however, see §7.3 for a discussion of the issues involved in 
calculating the asymmetric drift. 

5.3. Testing the outer exponential fit 

Finally, we test the hypothesis that the underlying disc scale length is given by the 
outer (> 2 kpc) portion of the SB profile by only fitting the interior points and the external 
points of the SB profile with a MCMC chain corresponding to scenario E. The outer cut is 
made at 2 kpc; Rd is subtly affected by the inclination and differences in scale height, and, 
moreover, the very tight correlation between Md and the disc Mj L makes it difficult to obtain 
reasonable acceptance rates with MCMC. Therefore, we manually tested combinations of Rd 
and M/L to produce the best fit; we adopt a value of Rd = 1.3 kpc and fix the disc M/L 
to 0.53Mrf. To avoid biasing the resulting value for Zd, we do not include the external disc 
in our \ 2 calculations. The inner cut is trickier; we have found that small changes to the 
number of inner points fit may lead to large changes in the resulting bulge parameters. To 
rectify this issue, we plot the residuals of the surface brightness from a regression line fit of 
Rd = 1.18 kpc (Fig. [TTI) and only fit the points that lie above the maximum of the Type 
II hump. At these points, the surface brightness is most assuredly bulge dominated. The 
stellar rotation is fit as in runs K and KL, along with the H I data in the outer part, and 
the velocity dispersion is fit as in the other three runs. 

For this run, we assume that the evolution of a disc would generate a ring of star 
formation that could reproduce the Type II hump. The process of generating this profile 
would alter the rotation and dispersions, however, and we caution that the kinematics that 
result may not fit the data as well as the direct fits presented here using MCMC. Accounting 
for this would require ad hoc selections of what points to fit, or excluding one or more data 
sets entirely; for example, excluding the stellar rotation and only fitting the H I data outside 
3 kpc. Such a process is probably no more reliable than simply fitting the full data sets, 
evolving a test model from the PDFs and then determining how well the rotation curve and 
dispersion profile fit the observed data after evolution. Because we are assuming that a bar 
is supposed to form in this scenario, unstable models would not rule out this interpretation. 
However, models with Q < 1 must be tested separately. 

As in scenario KL, the stellar rotation curve is well fit while the circular velocity is 
larger than the H/3 curve (Fig. [5]). The surface brightness fit is very good over the ranges 
we wish to fit; the inner eight points and outer exponential are well reproduced (Fig. [6]). 
The velocity dispersion profile looks different from that for the previous scenarios since the 
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bulge contributes significantly to the profile. Because of this, the a— drop is not properly 
reproduced (Sersic profiles do possess a slight central decline in the velocity dispersion - see 
Ciotti 1991 - but this feature cannot reproduce the observed a— drop). Thus, we posit that 
either an additional nuclear component is responsible for the a— drop, or that bar formation 
in this scenario would generate gas inflow that could reproduce the a— drop. The disc 
contribution to the velocity dispersion is comparatively muted, resulting in a lower overall 
disc dispersion. 

From Fig. [SJ we see that scenario E behaves differently than the other scenarios. These 
models have lower M d and Q m i n but Fig. [9] shows that these models have high A min . Addi- 
tionally, Q m in is reached at larger radii than in the other scenarios (typically, at ~ 2 — 2.5Rd), 
and they have more massive bulges (§6.4), which may also help inhibit the bar instability. 

As before, we selected several models and tested them for bar stability; the results are 
shown in the bottom right panel of Fig. [TUJ While these models quickly form spiral structure, 
bar formation is gradual at the low X min end of the best fit region and there is no evidence 
of bar formation at larger A min , despite the very low Q m i Q . Because these simulations do not 
include gas, it is not possible to simulate possible bar destruction by gas inflow, but we can 
examine the length of the bar that forms. In Fig. [12] we plot the bar strength A as a function 
of radius at different times for the lowest— X model. This plot provides an estimate of the 
bar length of ~ 1.5 kpc at 1.5 — 2 Gyr, which is consistent with the location of the Type II 
hump. We conclude that scenario E is a realistic scenario for the formation of the bar. 

6. Further results 

In this section we delineate numerous other results that can be gleaned from the model. 
Because scenarios KG and KL are ruled out, we focus mainly on the results for scenarios K 
and E. 



6.1. The disc 

We begin by comparing the photometric and dispersion scale lengths. We find from 
scenario K that the dispersion scale length is significantly larger than the photometric scale 
length, as seen in Fig. [131 The mean R a = 2.6 kpc is also considerably larger than the best 
fit R a of 1.16 kpc found in §2, while the photometric scale length is slightly less than the 'fit 
by eye' value of 1 kpc from B89. Note that the dispersion scale length is poorly constrained 
owing to the size of the error bars in the velocity dispersion data, and that a wide range 
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of R a may produce good fits, as Fig. [13] clearly shows. The large R a results from the low 
extrapolated central dispersion cq of only 37 km s -1 , much less than found by direct fit in 
§2. In scenario E, by contrast, the dispersion scale length of R G = 0.75 kpc is much lower 
than the adopted photometric scale length of Rd = 1.3 kpc. Here <7o is also very low, but 
the bulge dispersion is high enough to contribute much more to the dispersion profile than 
in the other scenarios, as is apparent from Fig. [7J Both scenarios suggest that <r and R a 
are not well constrained because of the large error bars in the LOS data. 

For the scale height z<i, we find a value of 0.14 ± 0.01 kpc in scenario K and 0.24 ± 0.02 
kpc in scenario E (Fig.ll4p. Both of these values are about one- fifth of the photometric scale 
length in scenarios K and E re spectively, and are consistent with observations of the scale 



height in edge-on disc galaxies (IKregel et al.l 120051 ) . 



The posterior PDFs for the hole radius Rh and the Kormendy index a are found in 
Fig. [15j We find Rh = 0.7 ± 0.1 kpc and a = 0.9 ± 0.1. Our value for a is in conflict 
with Kormendy's (1977) result that a = 3 for most galaxies (Baggett et al. 1998 also use 
a = 3 for their fits). A value of a = 3 yields a strong hole and therefore a very large bulge 
would be needed reproduce the SB profile. However, a large bulge would likely make the 
dispersion profile fit worse because the bulge is not likely to extend farther out than the 
highest dispersions at ~ 300 pc. We note that all three runs suggest a < 3; scenarios KL 
and KG yield a = 0.6 and a = 1.1, respectively. 



6.2. The bulge 

The results for the bulge parameters are found in Fig. [161 For scenario K, we find a 
bulge radius r& of 0.16 kpc and a Sersic index rib of 1.1, while scenario E yields a larger 
bulge with r& = 0.35 kpc and rib = 1.9. For scenarios KL and KG, the Sersic index is ~ 0.9. 
Scenarios K, KL and KG thus suggest that the bulge of NGC 6503 is nearly exponential. 
Note that the linear fit rib given in §2.1 is consistent with a pure exponential bulge. Scenario 
E sugg ests a cuspier bulge, but even here n b — 2, within the range of so-called 'pseudob- 



ulges' ([Kormendy fc Kennicuttl 120041 ) . Bulges of this type are thought to be generated via 
internal secular evolution processes - the gradual buildup of gas in the centre of the galaxy 
- rather than violent mergers (see, e.g., Kormendy et al. 2006). Given the NGC 6503 is 
isolated, secular evolution is likely the dominant evolutionary mechanism. These bulges typ- 
ically resemble small discs embedded in larger discs - their morphologies resemble flattened 
spheroids more than elliptical galaxies and their kinematics are rotationally dominated, but 
we cannot constrain the rotation of NGC 6503's bulge. 
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By way of comparison with known properties of bulges in spiral galaxies, the bulge 
radius is actually la rger than found by st udies of bulge-to-disc scale length correlation by 
roughly a factor of 2 ( iCourteau et al.lll996l ). However, NGC 6503's bulge also displays certain 
characteristics typical of nuclear cluste rs, even though its s ize exceeds typical nuclear clus- 
ters by about two orders of magnitude ( jWalcher et al.ll2005l ). Unlike normal bulges, nuclear 
clust ers are often offset from the dynamical centre of the galaxy OMatthews fc Gallagher 
20021 ). NGC 6503's bulge displays this phenomenon - the photometric cen tre is offset by 



100 pc. Nuclear bulges also have very small mass, typically <5x 10 7 M Q (jWalcher et al. 



20051 ); the mass of NGC 6503's bulge i s only slightly l arger . Thus, NGC 6503 has a surpris- 
ingly low density bulge. Furthermore, IWalcher et al.l ( 120051 ) find that their nuclear clusters' 
velocity dispersions are similar to that found in NGC 6503. Their nuclear M/Lj ratios are 
also consistent with the M/Lr values we find (§6.4), in that we expect the M/Li values to 
be slightly lower than our M/Lr values, which is precisely what occurs. Thus, NGC 6503's 
bulge properties are consistent with both nuclear clusters and ordinary pseudobulges. 

The o rigin of nuclear bulges is unclear, but it is plausible that they are formed by secular 
evolution. IWalcher et al.l (120061 ) find that repeated episodes of gas infall may contribute to 
nuclear cluster star formation, which is required to maintain their high luminosity. In the 
case of scenario E, bar formation may induce gas inflow that would generate a nuclear cluster 
that could in turn generate the a— drop. The most likely hypothesis is that of a pseudobulge 
that has a luminous nuclear cluster at the centre dominating the observed light. 



6.3. Halo cusps 

The halo cusp 7 is well-constrained, as Fig. [T7] attests. All runs suggest that the halo 
is cuspy, a result that is therefore quite robust. In the case of runs K and KL, we find 
7 = 1.2 and 7 = 0.9, respectively. For scenario KG, we find that 7 is close to 1 and 
therefore consistent with a NFW profile; this cusp is flatter than for the equivalent scenario 
that fits the stellar rotation curve. Meanwhile, scenario E suggest a very strongly cusped 
halo with 7 = 1.4. Thus, the derived cusp value depends on the density profile of the disc 
and the choice of rotation curve. Unsurprisingly, a stronger cusp is required to compensate 
for the lack of central mass in an inner truncated disc. In addition, we find a statistically 
significant difference in cusp values depending on which rotation curve is fit, but both cases 
are consistent with or cuspier than cosmological simulations. If the assumption that gas 
traces the circular velocity is erroneous, then the slope of the circular velocity curve must 
be steeper at the centre (as seen in Fig. [5]), and hence the derived cusp value will be larger. 
Thus, modelling asymmetric drift correctly is essential to obtaining the correct halo density 
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profile; we discuss possible issues with the asymmetric drift in §7.3. 



To verify that the halo is cuspy, we conducted a MCMC run in which the cusp value 
was fixed to 0. The resulting rotation curve fit is shown in Fig. HHJ the fit is clearly very 
poor. There is a pronounced kink in the rotation curve that is not observed, and neither the 
gas nor the stellar rotation curves are properly fit. We therefore reject models with 7 = 0. 

The issue of halo cuspiness remains unsettled. It is of t en argued that gala xies demon- 
strate evidence of cored halos (iBlais-Ouellette et al.l l200ll : Ide Blok et al.l 120011 ) , but many 
complications present themselves. Some researchers find that the rotation curves used to 
probe halo central structure are consis tent with both cored and cuspy profiles, in part be- 



cause of the size o f the uncertainties ( 



Dutton et al. 2005; Simon et al. 2005 



van den Bosch &: Swaters 2001: Swaters et al. 2003 



Spekkens et al.ll2005l ). Further caveats concern the 



correct interpretation of the rotation curve; the presence of noncircular gas motions, as might 
be caused by a triaxial halo, and the presence of turbulence caused by super nova feedback 



and the magnetorotational instability point to t he complexity of the p roblem (IHayashi et al. 



2004 : iRhee et al.l 12004 : IValenzuela et al.ll2007t iTamburro et al.ll2009l ). Some steps towards 
resolving these issues are taken by Oh et al. 2008, w ho find that the non-circular motions in 
their sample cannot rectify the discrepancy. In fact Begeman Jl987 ) finds t hat noncircular 
motions in H I are small in NGC 6503. However, ITamburro et al.l (120091 ) find relatively 
consistent values of ~ 5 — 25 km s _1 for the H I dispersion across their sample of late type 
galaxies, and the calculated star formation for NGC 6503 (§5.1) is consistent with their mea- 
surements. Therefore it is safe to assume that the H I dispersion of NGC 6503 is ~ 5 — 25 km 
s _1 and probably driven by supernova-induced turbulence inside the optical radius. More- 
over, as noted in §2.2, the H I rotation is systematically larger than the H/3 rotation, which 
suggests that the H/3 rotation is more affected by these issues. 

Recent simulations sug gest that halo cusps are slightly shallower than previously in- 
ferred, with 7 = 0.9 ± 0.1 ( INavarro et al.l 120081 ). However, halos are subject to multiple 
effects due to baryons that are not accounted for in cosmological simulations, which may 
help explain why our cusp values are larger than found by Na yarro et al. The fu ll effect of 
baryonic physics on halo profiles is still poorly understood, but lAbadi et al.l (120091 ) find that 
the central density increases when more detailed b aryonic effects are includ ed. The central 
density cusp may also be affected by bar formation; [Weinberg fc Katzl (120021) argue th a t bars 
are capable of washing out c entral cusps, an effect seen in iHolley-Bockelmann et al.l (120051 ) 
and IWeinberg fc Katzl ( 120071 ). However, Sellwood (2003, 2008) argues that bar formation 
can draw mass inward, increasing the strength of the cusp; iDubinski et al.l (120081 ) find with 
their simulations that a live bar maintains the cusp. The differences in cusp evolution may 
be due to differences in the codes used. If NGC 6503 did once possess a bar, our results 
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argue against a bar-induced cusp flattening in this galaxy. 



6.4. Masses and M/L ratios 



The halo mass we find is necessarily limited by the outermost H I data point at 800", and 
so should be regarded as a lower limit. The outer rotation curve is the primary determinant 
of the halo mass, so we do not expect significant variations in halo masses across the different 
runs. We find M 2 o ~ 60 x 1O 9 M across all scenarios. For scenario K we obtain a = 
3.2 ± 0.2 x 10 9 M Q and M b = 7.0 ± 1.4 x 1O 7 M . Larger disc masses are obtained for 
scenario KL, while the bulge mass is much lower for scenario KG. In scenario E, we find 
M d = 3.0 ± 0.4 x 10 9 M o and M b = 15 ± 2 x 1O 7 M . Unsurprisingly, in all cases the halo 
accounts for > 90% of all the galaxy mass, while the bulge mass is too small to significantly 
impact the stability of the disc, except in scenario E. Whether or not the bulge masses 
are consistent with other values found in the literature depends on whether the bulge is 
interpreted nuclear cluster (§6.2). 



Mj L ratios for the disc and bulge are shown in Figs. [19] and [20j In order to reproduce 
the u — drop, the bulge M/L ratios must be small, as we find in scenarios K, KL and KG, 
but not E. The first three scenarios show that the bulge M/L is less than 1, while the disc 
M/L is roughly twice that (for scenario KG, it is an order of magnitude larger). We find the 
disc M/Lr to be 1.9 ± 0.3 and the bulge M/Lr to be 0.9 ± 0.1 under scenario K. In scenario 
E, we fixed the disc M/Lr to be 0.53M^, yielding a disc M/Lr of 1.6 ± 0.2 while the bulge 
M/Lr is 1.2 ± 0.1, which, although larger than 1, is still lower than the disc M/Lr. 

Our models assume constant Mj L throughout the galaxy disc in the absence of a stellar 
population gradient. In parti cular, stellar M/L ratios tend to increase in s piral galaxy centres 
where redder colours prevail ( IBell &: de Jong||2001t Ide Jong fc Belli 120071 ). From B89, B — R 
across two disc scale lengths for NGC 6503 varies from 1.10 t o 1.41, correlating with M/Lr 
ratios of ~ 1.3 — 2.7 using Table 3 of IBell fc de Jongj (120011 ) . which is consistent with our 
findin gs. Upper limits to the dynamical M/ L R values were obtained by iBroeils fc Courteau 
( 119971 ). who find M/L R < 4, also consistent with our values. 



The M/L values for the bulge do not, by c ontrast, coincide with othe r bulge values 
found in the literature for nearby spiral galaxies. lYoshino fc Ichikawal ( 120081 ) generally find 
that their bulge M/L ratios are larger than their disc values with few exceptions; M/Ly 
ratios are ~ 4.5 ± 2.4 and M/Lj ratios are ~ 2.7 ± 1.8. In no case is the bulge M/Ly value 
lower than 1, and only in a few cases does M/Lj fall below 1. NGC 6503's bulge M/L lies 
significantly below nearly all the bulges in this sample. However, M/L ratios are highly 
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sensitive to the fo rmation of mas s ive, l umino us stars, and our small M/L v alues, together 
with the work of IWozniak et al.l (120031 ) and IWozniak fc Champavertl (120061 ) showing that 
o — drops can be caused by massive star formation, strongly suggest that the bulge of NGC 
6503 has a star- forming c ompo nent. Thus, the better comparison may be to the nuclear 
clusters of IWalcher et al.l ( 120051 ). who find M/Lj ratios that generally lie below 1. 



7. Discussion 

7.1. Comparison with earlier work 

B89 and BG97 present dynamical models for NGC 6503, and BG97 test the models for 
bar stability. In both papers, the mass distribution of the disc is given by Eq. CGl] with a 
scalelength of Rd = 1 kpc, while the velocity ellipsoid is given by 



a z = ^7rE(R)z d (19) 

vr = Vz/0.6 (20) 

^ = a Ry /B/(B - A) (21) 

where E(i2) is the project ed surface density, Zd is the scale height of the disc, and A and 



B are the Oort constants (Ivan der Kruit fc Searldll98ll ). The expression for a z is obtained 



for an isolated isothermal disc, while the relationship between a R and comes from the 
epicycle equations. However, the relationship between <jr and a z is an extrapolation from 
the observed value in our solar neighbourhood. We return to this assumption in §7.2. 

B89 also examined a s econd model in which Q is held constant throughout the disc 



(jCarlberg fc Sellwoodlll985l ). For this model, the radial dispersion obeys 

<t r <x £ exp(-R/R d ) (B(B - A))' 1 ' 2 . (22) 



K 



Both models yield reasonable fits (the former is slightly better ), but cannot reproduce th e 



u — drop. To reproduce a borderline stable disc where Q = 1.7 ( ISellwood fe Carlberglll9841 ) 



B89 finds that the disc M/L must be 1.7 ± 0.3 in the B— band, which is consistent with our 
findings. 

BG97 model the galaxy in more detail and provide stability studies. They adopt the 
gas rotation curve as the fundamental input to their simulations and assume an isothermal, 
cored halo profile along with the disc from B89. They vary the disc to halo mass ratio to 
investigate the stability properties of the galaxy; Q declines slightly from the lowest disc to 
halo mass ratio to the highest. 
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The dispersions of the isolated disc scale with the square root of its surface density 
(Eqs. [19] and [20]) , but the process of embedding the disc in a dark halo modifies the disper- 
sions. BG97 find that embedding the isolated disc in the halo allowed the settling process at 
the beginning of the N— body simulations to rearrange the dispersions itselfjfl The result is 
an initial central outflow in the simulations, leading to a decrease in the dispersions, more- 
or-less naturally accounting for the observed disc dispersion of NGC 6503. BG97 find that 
their two lightest discs remain stable to bars, but are unable to reproduce the a— drop. 

Visible on Fig. [S] are the points corresponding to the models used by BG97. They are 
located well above the 2cr confidence interval. The discrepancy occurs because Q is obtained 
from the radial dispersion ctr, which, in BG97's models, is directly coupled to the surface 
density (cf. Eqs. [19] and [20]) . Thus, BG97's ctr is determined by the disc surface density, 
which is not the case for the GalactlCS model. Their Fig. 9 shows that, after settling, 
their minimum Q— values decline so that they fall in line with our PDF in Fig. [8] Thus, 
their Md — Qmm data points are consistent with our Md — Q m m PDF. However, their discs 
do not display bar formation, while ours do; there are two possible reasons for this. First, 
from Fig. [TD], most of our models only develop a bar after ~ 2 Gyr, while BG97 evolved 
their simulations for ~ 1.3 Gyr. Second, and more i mportantly, it is now known that live 
halos can trigger bar formation (lAthanassoulal 120021 ) ; BG97 used rigid halo potentials for 
their simulations, so this bar formation trigger is not in play. We evolved one of our bar 
unstable models in exactly the way described in §5, except the number of disc particles was 
reduced to 40K. We found that the bar formation in this model was delayed by ~ 1 — 2 Gyr, 
suggesting that sufficient numerical resolution is needed to capture the correct evolution. 
This is consistent with halo-triggered bar formation, because inadequate numerical resolution 
will fail to resolve the resonant interactions required for this mechanism. 



3 BG97 attempted to account for the embedding: a z can be modified by a factor F(e) accounting for the 
relative halo contribution, where 

The result is that when e is less than (as it could be near the centre where -§rv1 is very large) the 
dispersion can decrease upon embedding, while large values of e cause the dispersions to increase. BG97 
dropped this factor for the models discussed in that paper and used the isolated disc dispersions (Bottema, 
private communication) . 
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7.2. The ratio of velocity dispersions 



The ratio <y z /a R = 0.6 is an observatio n of the solar n eighb ourho od. Should th i s ratio 
hold at all radii in a general disc galaxy? iGerssen et al.l (119971 ) and IGerssen et al.l (12 000) 



find that the ratio is closer to 0.7 for NGC 488 and to 0.85 for NGC 2985, and Westfall et al. 



()2007j) find that the ratio is also high for NGC 3949 and NGC 3982 (the latter's vertical 
dispersion is larger than its ra d ial dis pers ion) . Theoretical studie s of the dispersion ratio have 
been carried out by llda et al.l (119931 ) and IShiidsuka fc Ida (jl999l ). who find that (TR/cr z ~ 0.6 
is roughly correct provided that k/Q < 1.5; higher values of k/VL imply higher a z values. 
This condition is probably satisfied in the centres of galaxies, and it would therefore be 
unreasonable to extend a fixed ratio of velocity dispersions across the entire disc. In fact, 
one might surmise that the reason for the central outflow in BG97's heavier discs is the 
incorrect ratio of velocity dispersions there. 

To examine the velocity dispersion ratio for NGC 6503, we took the fiducial best fit 
model for NGC 6503 (that is, the means from Table W) and all models used to test bar 
stability and calculated the ratio as a function of radius. The results are found in Fig. [2TJ 
NGC 6503's dispersion ratio depends critically on the assumed model - scenario K exhibits a 
nearly linear decline with radius in which a z /<Jn > 0.6 within two disc scale lengths. There is 
no indication of flattening that would justify adopting a single value over the entire disc, and 
very little scatter. Scenarios E, KL and KG display considerably more scatter, although only 
scenario KG displays anything close to a relatively constant ratio of 0.6 over a significant 
range; models from scenario KG appear to converge to o z ja R ~ 0.6 at small radii, but the 
scatter increases at larger radii. The shape and slope of the curve is critically dependent on 
the choice of Rd and R a ; in scenario K, R d <C R ai which drags down the dispersion ratio 
at outer radii. If we set R a to be 50% of the fiducial value, the ratio would fall less rapidly 
with radius but would still not be constant. Note that when Rd = R a the ratio is constant 
by Eqs. [19] and |20j although this may not hold in an embedded disc. Overall, we find that 
a constant dispersion ratio does not apply for this galaxy. 



7.3. The relationship between asymmetric drift, stability and cusp value 

The issues of cusps, asymmetric drift and stability are intertwined because all hinge on 
the reliability of rotation curve data and all depend on a full accounting of effects that alter 
the interpretation of the rotation curve. As already indicated, it is inappropriate to assume 
that the gas rotation traces the gravitational potential, implying that: (1) the difference 
between the gas and stellar rotation curves is not a measure of the asymmetric drift; (2) the 
gas rotation cannot be used to assess the cusp value; and (3) the proper choice of rotation 
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curve to model is critical to finding the correct stability region. As noted in §2.2, the 
theoretical asymmetric drift (Eq. |3]) is much larger than the difference between the gas and 
stellar rotation curves. We have conducted separate MCMC runs in which we fit both the 
gas rotation curve as a tracer of the circular velocity and the stellar velocity. This only yields 
models with Q < 1, because of the very low asymmetric drift implied by the nearly coincident 
rotation curves. Such models, when simulated, display sharp instabilities including disc 
fragmentation and cannot reproduce the observed properties of the galaxy. Our preference 
for fitting the stellar rotation curv e and finding t hat g as cannot be assumed to trace to 
circular velocity is corroborated by iPizzella et al.l ( 120081 ). They analyse rotation curves of 
LSB galaxies and find that stellar rotation curves are much more regular and amenable 
to modelling than gas curves, which suffer from numerous issues including noncircular and 
vertical motions that affect their speeds relative to the circular velocity. 

The impact on the cusp value is notable because, if gas does not trace the gravitational 
potential, the circular velocity is steeper in the centre than the slope of the gas rotation. 
Assuming that the gas traces the gravitational potential will therefore lower the inferred 
cusp value. This is a possible source of error in early work assessing the observed cusp value 
in LSB galaxies, as small changes in central slope can lead to large changes in the cusp 
value or possibly imply a cored halo (7 = 0). We avoid this problem by using a model that 
follows the stellar kinematics with asymmetric drift incorporated self-consistently. Once the 
gas asymmetric drift is properly accounted for by including a full treatment of noncircular 
motions and turbulence, gas rotation curves should be as reliable as stellar rotation curves. 
Detailed measurements of the H I dispersion in spiral gal axies are becoming mo re common 
and will aid the modelling of gas kinematics. For example, iBoomsma et al.l (120081 ) obtain the 
high resolution H I velocity map for NGC 6946, finding dispersions of ~ 6 — 13 km s _1 , high 
velocity H I clumps that lag the disc rotation, and hundreds of "holes" in the H I distribution 
which are likely due to star formation. We surmise that modelling these phenomena and 
determining the impact on derived galaxy parameters is a nontrivial exercise. However, while 
easier to model than the gas asymmetric drift, the stellar asymmetric drift is also subject to 
several complications, a few of which we now address. 



Eq. [3] is derived assuming an exponential disc, cylindrical alignment of the velocity 
ellipsoid and Rd = R a . With an inner truncated disc and Rd 7^ R a , the asymmetric drift 
becomes 



v a 



2Vr 



R R 1 (Rdv s 
R~d + R~ a + 2\V s dR 



1 



1 
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Rh 
R 



(24) 



If we now assume that the velocity ellipsoid is spherically aligned, we obtain 



a 



R 

2v, 



R R 1 / R dv s 
R~ d + R^ + 2\V s dR 
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R 



(25) 
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This expression will match Eq. |23]only if & z = ctr, which is not generally the case (§7.2). 
Both expressions reduce the asymmetric drift relative to what Eq. [3] predicts; in particular, 
the addition of an inner truncation reduces the expected drift significantly in the inner parts. 
The size of the change is only a few km s _1 , but it is still a significant fraction of the predicted 
drift. This is seen in Fig. [221 i n which the effects of adding an inner truncation and spherical 
alignment of the velocity ellipsoid are examined. It is noteworthy that the inner truncation 
can effectively wipe out the asymmetric drift over a large range in radius; a similar effect is 
seen in Fig. [5j For Fig. [22} we use the fiducial values for scenario K found in Table El Note 
that the increased asymmetric drifts at larger radii are due to the large R a of this scenario, 
and the asymmetric drift expression breaks down inside ~ 0.5Rd- We test cases in which R a 
is reduced to 50% of the fiducial value and find that the asymmetric drift is reduced in the 
outer part. Thus, as is evident from Fig. [22l a combination of factors can serve to reduce the 
asymmetric drift nearly to the observed values in this galaxy. One should keep in mind that 
the sensitivity of the asymmetric drift expression to such changes means that the predicted 
drifts may fluctuate wildly depending on the assumptions. 

Because of these issues, the circular velocity derived in scenarios K, KL and E may be 
larger than the true circular velocity. These caveats further suggest that the asymmetric 
drift found in scenario KG may be smaller than seen in Fig. [51 and hence that scenario 
KG may be potentially viable. Without more information about the asymmetric drift in 
real galaxies, it is impossible to claim that scenario KG is definitely unrealistic; however, 
scenarios K and E remain our preferred models for this galaxy. 



7.4. Evidence for a bar 



Recent work by iFreeland et al.l (120101 ) suggests that NGC 6503 may possess a nuclear 
disc of radius ~ 100 pc and an end-on bar; deprojection of their if— band image of NGC 
6503 shows a bar like structure perpendicular to the major axis with clear spiral arms 
emanating from the ends. This scenario is consistent with our scenario E. It would moreover 
explain the a— drop because o f star formation at the centre ca used by bar induced gas inflow 
and because diagnostics from iBureau fc Athanassoulal (120051 ) suggest that end-on bars may 
generate a— drops. These findings are compelling evidence for a bar in NGC 6503. Notably, 
these authors also estimate the age of the star format ion ring at ~ 0. 5 Gyr , which is well 
below the timescale for bar destruction determined by iBournaud et al.l (120051 ) of ~ 1 Gyr. 



The presence of a bar in NGC 6503 is unlikely to change our main conclusions because 
it is end on. There is no clear kinematic indicator of bar structure in the stellar or ionized 
gas rotation curve; we expect that the measured values are largely unaffected by the bar 
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since only a small region of the slit used for observation would have intersected the bar. Part 
of the & — drop could be due to the bar, and thus the bulge may not be as large and have 
as low a density as found in §6.2. Other bulge parameters are less likely to be affected by 
the bar; in particular, the conclusion that there is a star forming component in the galaxy 
centre is unchanged. 

8. Conclusion 

We have deployed a Bayesian/Markov chain Monte Carlo technique to model the disc 
galaxy NGC 6503. We find models for NGC 6503 that satisfy the observational constraints 
using this technique, and are able to constrain the input parameters. The data includes a 
Freeman Type II surface brightness profile, ionized gas and H I rotation curves, the stellar 
rotation curve and the stellar line of sight velocity dispersion. Four different scenarios were 
considered: an inner truncated disc, an exponential disc with dust, a model in which gas 
traces the gravitational potential, and a model in which the true underlying disc scale length 
is revealed in the outermost portion of the surface brightness profile. We find that the second 
scenario leads to models that strongly bar unstable, and the third scenario cannot reproduce 
the stellar rotation curve. The first scenario provides the best fit to the data and is most 
likely to offer the bar stability required to correctly model the galaxy, but the last scenario 
also provides a realistic model for the galaxy. 

Further properties of the galaxy are discerned; the bulge is a pseudobulge, and the bulge 
M/L is lower than the disc M/L, suggesting a star forming component that is probably 
responsible for the a— drop. We also find that the halo must be cusped with 7 > 1, a result 
that is robust to all fitting methods. The Bayesian/MCMC technique used to discover these 
results is robust and flexible; we find it to be an effective tool for fitting galaxy models in 
complex parameter spaces. 

We thank K. Spekkens, P. Teuben, R. Bottema, J. Dubinski, M. Bershady and C. 
Arsenault for useful and informative discussions. S.C and L.M.W. wish to acknowledge 
support through respective Discovery Grants from the Natural Sciences and Engineering 
Research Council of Canada. 
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Fig. 1. — Surface brightness and colour profile from Bottema (1989). The black curve shows 
a linear fit to the outer exponential with scale length 1.18 ± 0.04 kpc. 
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Fig. 2. — A log/i-logi? plot of the surface brightness. The green line is determined by the 
slope of the inner two points and has a slope of ~ 1. The black curve shows the fit found in 
Fig. CD The central magnitude fiQ is estimated from the slope of the inner points in Fig. CD to 
be 17.6 mag arcsec -2 . 
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Fig. 3. — Rotation curve data for NGC 6503. The H/3 data are filled red circles, the stellar 
rotation data are filled blue triangles (Bottema 1989) and the H I data are hollow purple 
squares (Begeman 1987). The black curve shows the fit using Eq. [2] to the gas data and 
the green curve shows the fit to the stellar data. The fitting parameters are: for the black 
curve, vq = 117 ± 1 km s _1 , R c = 0.88 ± 0.03 kpc and £ = 2.1 ± 0.1; and for the green curve, 
v Q = 108 ± 4 km s~\ R c = 0.95 ± 0.04 kpc and ( = 3.2 ± 0.6. The data have not been 
corrected for inclination. 
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Fig. 4. — LOS stellar dispersion data from Bottema (1989) in red. The black curve shows 
an exponential + inverse Gaussian (Eq. H|) fit to the data. The fitting parameters are given 
by a = 59.5 ± 8.4 km s~\ R a = 1.16 ± 0.29 kpc, A = 0.59 ± 0.06 and B = 0.12 ± 0.02 kpc. 
The data have not been corrected for inclination. 
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Fig. 5. — Graphs of the rotation curves for representative models of each scenario. Red 
identifies the H/3 points, purple identifies the H I data and blue represents the stellar rotation 
points as in Fig. |3j Green is the stellar model circular velocity curve while black represents the 
circular velocity curve. The orange curves show the rotation curve breakdown by component: 
halo (solid), disc (dashed) and bulge (dotted). The scenario K and KG rotation curves display 
a kink in the black curve at ~ 0.2 kpc, owing to the inner truncation. 
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Fig. 6. — Graphs of the SB profile for scenarios K and E. Orange is the bulge contribution 
to the surface brightness, green is the disc contribution, black is the total and red identifies 
the observed R— band surface brightness. All models from scenarios K, KL and KG show 
fits similar to the top panel. 
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Fig. 7. — Graphs of the LOS dispersion for scenarios K, KG, and E. The model dispersion 
is in black and the observed dispersion is in red. The bulge contribution to the dispersion is 
in orange and the disc contribution is in green. Fits for scenario KL are similar to those for 
scenario KG. 
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Fig. 8. — The 2-D PDF of the disc mass versus the minimum Q, Q m i n , of the disc. Black 
refers to scenario K; green refers to scenario KL; red refers to scenario KG; blue refers to 
scenario E. The thick lines enclose the la confidence region, and the thin lines enclose the 
2a confidence region. The blue-outlined stars correspond to the parameters of the models 
chosen to test stability, and the purple circles correspond to the points used by BG97 (§7.1). 
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Fig. 9. — The 2-D PDF of the minimum X, X min , versus minimum Q, Q m i n , of the 
Colours and symbols are as in Fig. [8j 
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Fig. 10. — The Fourier amplitude In A of several models obtained from each scenario. The 
colour pattern orange-green-blue-red corresponds to lowest to highest Q ui \ n of the models 
identified by stars on Figs. |8] and [9j 
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Fig. 11. — Residuals in the R— band surface brightness for the exponential disc fit shown in 
Figure [TJ 
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Fig. 12. — Bar strength A as a function of radius at different times for the lowest X model 
from scenario E. Pink is the bar strength at 1.5 Gyr, purple at 2 Gyr, orange at 3 Gyr, and 
green at 4 Gyr. 
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Fig. 13. — The 2-D PDF of the photometric scale length Rd versus the dispersion scale 
length R a . Colours are as in Fig. [SJ The dotted line identifies the Rd = R a line. 
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Fig. 14. — The PDF for the disc scale height Zd- Black refers to scenario K; green refers to 
scenario KL; red refers to scenario KG; blue refers to scenario E. 
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Fig. 15. — The 2-D PDF of the Kormendy cutoff radius Rh versus the cutoff index a. Colours 
are as in Fig. [SJ 
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Fig. 16. — The 2-D PDF of the Sersic index n& versus the radius of the bulge r&. Colours 
are as in Fig. [SJ 
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Fig. 17. — The PDF for the halo cusp strength 7. Colours are as in Fig. [TU 
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Fig. 18. — Rotation curve for a model from a MCMC run in which the cusp is fixed to 0. 
Colours are as in Fig. |5j 
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Fig. 19. — The 2-D PDF of the disc mass M d versus the disc mass-to-light ratio (M/L) d . 
Colours are as in Fig. |HJ 
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Fig. 20. — The 2-D PDF of the bulge mass M b versus the bulge mass-to-light ratio (M/L) b . 
Colours are as in Fig. |HJ 
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Fig. 21. — Ratio of velocity dispersions for NGC 6503. Black refers to scenario K, green 
refers to scenario KL, red refers to scenario KG and blue refers to scenario E. The thick lines 
use the best fit values found in Table [2] for scenario K and the best fit values for scenario KG 
(not published), while the thin lines correspond to the specific models identified by stars on 
Fig. |8j There is no thick line for scenario KL because the means for scenario KL typically 
fall between two separate PDFs (§5). 
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Fig. 22. — Calculated asymmetric drifts in the case of cylindrical alignment and Kormendy 
hole (solid blue), cylindrical alignment and exponential disc (dashed blue), cylindrical align- 
ment, Kormendy hole and low R a (dotted blue), spherical alignment and Kormendy hole 
(solid orange), spherical alignment and exponential disc (dashed orange), and spherical align- 
ment, Kormendy hole and low R a (dotted orange), using the fiducial parameters for scenario 
K found in Table |5J The black points are the observed asymmetric drift, given by the 
difference between the H/3 rotation and the stellar rotation. 
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Parameter 


Units 


Description 


r h 


kpc 


Halo truncation radius 


Vh 


100 km s- 1 


Halo characteristic velocity- 


a-h 


kpc 


Halo characteristic radius 


5 


kpc 


Halo truncation width 


7 


Dimensionless 


Cusp value 


M e 


1O 9 M 


Exponential disc mass 


Rd 


kpc 


Photometric disc scale length 


z d 


kpc 


Disc scale height 


n b 


Dimensionless 


Sersic index 


v h 


100 km s- 1 


Bulge characteristic velocity 


n 


kpc 


Bulge characteristic radius 




100 km s- 1 


Central velocity dispersion 


^rot 


Dimensionless 


Bulge rotation parameter 


(M/L) d 


M & /L Q 


Disc mass to light ratio 


(M/L) b 


Mq/Lq 


Bulge mass to light ratio 


R c 


kpc 


Radial dispersion disc scale length 


Rh 


kpc 


Kormendy cutoff radius 


a 


Dimensionless 


Kormendy cutoff index 


Calculated Quantities 


M d 


1O 9 M 


Real disc mass 


M b 


1O 9 M 


Model bulge mass 


M 20 


1O 9 M 


Model halo mass within 20 kpc 


Q 


Dimensionless 


Toomre local stability parameter 


X 


Dimensionless 


Toomre global stability parameter 



Table 1: Table of the input parameters for the MCMC chains, as well as calculated quantities 
that are output for each model. 
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Parameter Best-fit value 



Observed value 





Scenario K 




Scenario E 








2.47 


±0.02 


2.38 


±0.02 




a h 


6.66 


±0.49 


8.69 


±0.66 




7 


1.22 


±0.04 


1.43 


±0.03 




M e 


5.77 


±0.83 


3.12 


±0.43 




Rd 


0.71 


±0.02 


1.3 


fixed 


1 




0.14 


±0.01 


0.24 


±0.02 


< 0.2R d 


n b 


1.06 


±0.07 


1.93 


±0.10 




Vb 


0.43 


±0.04 


0.51 


±0.03 




T i 


0.16 


±0.02 


0.35 


±0.02 






0.37 


±0.03 


0.34 


±0.03 


0.55 ±0.10 


{MIL) 


d 1-86 


±0.33 


1.60 


±0.2 


1.41 ±factor of 2 


{MIL) 


6 0.87 


±0.12 


1.23 


±0.11 




R a 


2.57 


±0.29 


0.75 


±0.05 


1 


Rh 


0.69 


±0.11 








a 


0.90 


±0.14 








M d 


3.19 


±0.59 


3.03 


±0.39 




M b 


0.070 


±0.014 


0.15 


±0.02 




M 20 


61 


±2 


61 


±2 




Q 


1.62 


±0.24 


0.69 


±0.15 


1.3 


X 


1.88 


±0.46 


2.84 


±0.57 




Table 2: 


The final posterior 


values for each input and calculated parameter with la error 



bars, with units as in Table [TJ Where possible, the observed values are obtained from B89 
or (in the case of the disc mass to light ratio) BG97. 
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